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Abstract 

\ 

^This  thesis  presents  the  design  of  a proposed  terrain 
following  flight  controller  capable  of  tracking  a cubic 
spline  reference  path.  The  controller  tracks  only  the 
longitudinal  motion;  however,  it  provides  both  thrust  and 
elevator  control. 

The  controller  is  based  on  a nonlinear  reference 
model  constructed  from  an  optimized  spline  path.  The 
spline  path  altitude  and  associated  derivatives  are 
computed.  The  longitudinal  equations  of  motion  are  linear 
ized  about  the  reference  trajectory  and  a truth  model  of 
actual  aircraft  motion  is  developed.  Deviations  from  the 
nonlinear  reference  trajectory  to  the  actual  provide  a 
linear  system.  Optimal  control  theory  is  used  to  solve 
the  resulting  linear  regulator  problem.  The  feedback 
gains  are  calculated  from  the  steady-state  Ricatti  matrix 
equation.  The  two  system  controls  are  updated  using  these 
feedback  gains. 

Results  show  the  feedback  controller  to  be  stable 
and  capable  of  tracking  the  cubic  spline  reference  path. 
The  performance  index  weighting  matrices  can  be  adjusted 
to  improve  the  controller.  The  resulting  controller  can 
provide  both  good  path  control  and  engine  control  to 
improve  vehicle  survivability,  engine  life  and  fuel 
consumption. 


OPTIMAL  TERRAIN  FOLLOWING  CONTROLLER  FOR 


AN  OPTIMIZED  SPLINE  REFERENCE  PATH 

I . Introduction 

Lov;  altitude  high-speed  flight  has  become  a necessity 
to  increase  penetration  probability  and  survivability  of 
strategic  bomber  aircraft.  This  has  resulted  in  emphasis 
being  placed  on  the  terrain  follov/ing  mission  profile. 

And  in  recent  years,  the  terrain  following  field  has 
broadened  to  include  strategic  cruise  missile  systems  and 
tactical  helicopters.  This  broadening  has  been  accompanied 
by  the  incorporation  of  all  the  best  features  of  modern 
technology. 

Current  Problem 

The  terrain  follov;ing  mission  creates  severe  demands 
on  the  propulsion  system  of  an  aircraft.  These  demands 
take  the  form  of  large  and  frequent  fluctuations  in  the 
thrust  control  system  that  result  in  shorter  engine  life 
and  increased  fuel  consumption.  The  operation  of  the 
F-111  aircraft  during  the  last  ten  years  bears  out  this 
problem.  Maintenance  reports  conclude  that  the  terrain 
following  mission  adversely  affects  both  engine  wear  and 
fuel  consumption  (Ref  2:4). 

Modern  control  theory  has  progressed  to  the  point 
where  effective  computational  techniques  can  definitely 
minimize  the  disadvantages  of  terrain  following  flight. 
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Also,  the  state  of  the  art  in  on-board  computer  systems 
is  such  that  the  required  real-time  calculations  are 
achievable.  The  coupling  of  control  theory  with  better 
computation  leads  to  development  of  better  terrain  follow- 
ing systems  which  reduce  engine  wear  by  minimizing  throttle 
control  motion. 

Currently,  the  problem  which  exists  is  the  design  of 
a system  that  optimally  controls  all  the  aspects  of  ter- 
rain following  flight.  The  controller  must  be  capable 
of  maintaining  minimum  deviations  from  a designated 
reference  path  while  providing  reduced  engine  wear  and 
decreased  fuel  consumption. 

Background 

The  fundamental  problem  in  terrain  following  flight 
is  determining  the  commands  to  the  control  system  which 
achieve  a desired  flight  path.  Extensive  research  has 
been  addressed  to  this  one  area  alone.  As  a result,  many 
methods  and  computer  algorithms  have  been  developed  to 
fly  the  terrain  following  mission.  These  methods  encom- 
pass the  flight-path  angle  controllers  and  the  path 
controllers.  Each  of  these  categories  contains  a hierarchy 
of  improved  development  featuring  an  optimal  tracking 
system  at  the  top.  Accounting  for  the  present  state  of 
the  art,  these  optimal  trackers  provide  accurate  naviga- 
tion with  minimal  deviation  from  the  nominal  flight  path. 
But  none  of  these  optimal  methods  adequately  accounts  for 
the  thrust  control  problem. 
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The  t.'rust  contro]  problem  has  not  been  considered 
because  present  terrain  following  control  laws  essentially 
decouple  vertical  flight  path  control  from  speed  control. 
These  control  laws  assume  a constant  speed  as  the  aircraft 
contours  the  terrain.  It  can  be  seen  that  maintaining  a 
constant  speed  throughout  a wide  variety  of  flight-path 
angles  will  certainly  cause  large  engine  thrust  fluctua- 
tions. Even  though  present  terrain  follov/ing  control  laws 
use  the  constant  speed  assumption,  recent  tests  have 
shown  that  the  overall  system  still  performs  well  when 
the  actual  speed  varies  slightly  (Ref  8:1).  Apparently, 
the  control  system  is  fairly  insensitive  to  reasonable 
speed  variations.  Therefore,  a control  law  v/hich  allows 
reasonable  velocity  deviations  from  nominal  will  produce 
lessened  throttle  adjustment.  The  best  way  to  achieve 
this  is  to  incorporate  throttle  control,  speed  of  the 
aircraft,  and  vertical  flight  path  so  that  these  all  will 
interrelate  and  be  optimally  controlled.  This  allows  a 
trade-off  betv/een  speed  variation  and  thrust  fluctuation 
while  flying  the  desired  reference  flight  path. 

The  Boeing  Aerospace  Company  attacked  the  problem 
using  an  optimal  integrated  terrain  following  system 
(Ref  8).  This  optimal  system  controls  both  altitude  and 
speed  while  accounting  for  vertical  acceleration  con- 
straints and  available  engine  thrust.  Boeing’s  design 
minimises  a quadratic  performance  index  which  includes 
the  following  characteristics: 


'j. 
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(a)  Altitude  deviation  from  desired  flight  path 

(b)  Speed  deviation  from  set  speed 

(c)  Vertical  acceleration  command 

(d)  Throttle  command 

The  adjoint  differential  equations  for  the  feed  forward 
tracking  signals  are  processed  in  accelerated  reverse 
time  to  generate  the  current  terrain  following  commands. 
Once  each  second  a now  desired  path  is  computed. 

The  results  of  the  Boeing  integrated  model  show  that 
the  throttle  motion  was  reduced  dramatically  when  compared 
to  current  terrain  following  control  systems.  Only  half 
the  excess  fuel  used  by  previous  controllers  was  consumed 
on  the  same  profile.  The  Boeing  simulation  maintained 
aircraft  speed  closer  to  the  set  speed  and  minimised  the 
altitude  standard  deviation  over  a typical  flight  profile. 
The  only  adverse  effect  was  that  the  aircraft  vertical 
acceleration  was  increased,  and  this  resulted  in  a degra- 
dation of  ride  quality.  This  is  quite  important  for  pilot 
comfort,  and  can  be  imp»'oved  by  adjusting  the  acceleration 
limits  within  the  system  model  (Ref  8:46), 

Overall,  the  Boeing  design  provides  the  basis  for 
an  optimal  controller,  but  it  does  not  provide  for  an 
op.timizod  reference  path.  The  optimised  cubic  spline 
reference  path  generator  developed  by  Major  James  E.  Funk 
achieves  this  goal.  The  optimised  spline  path  incorpo- 
rates clearance,  slope  and  acceleration  constraints  in 
the  optimisation  problem,  and  generates  a cubic  spline  as 
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a reference  path.  This  path  has  continuous  first  and 
second  derivatives.  Each  solution  to  this  optimization 
problem  results  in  a smooth  reference  path  v/hich  follows 
the  terrain  as  closely  as  possible.  Because  this  smooth 
path  contains  continuous  first  and  second  derivatives, 
an  aircraft  should  easily  be  able  to  fly  this  path.  So, 
if  an  optimal  controller  can  be  found  which  uses  the 
optimised  cubic  spline  generated  path  as  a reference,  the 
aircraft  controllability  will  be  maximized  and  the  thrust 
control  problem  will  be  minimized. 

Scope 

The  object  of  this  thesis  is  to  demonstrate  the 
feasibility  of  an  optimal  controller  which  attempts  to 
fly  an  optimized  spline  reference  path.  This  report  uses 
optimal  control  theory  v;ith  special  emphasis  on  optimal 
linear  state  feedback  control. 

The  optimal  controller  is  based  on  deviations  from 
a nonlinear  reference  model.  The  feedback  gains  are 
computed  for  a system  linearized  about  the  reference 
model  based  on  a quadratic  performance  index  somewhat 
different  from  the  one  used  by  Boeing.  The  linearized 
equations  are  also  used  to  estimate  the  proper  value  for 
the  weighting  matrices.  The  feedback  gains  will  be  used 
in  the  nonlinear  simulation  of  aircraft  motion  to  evalu- 
ate the  optimal  controller.  Performance  of  the  controller 
v/ill  be  evaluated  under  varying  conditions. 

This  study  v;ill  investigate  the  optimal  controller 
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for  aircraft  only.  The  aircraft  data  used  in  the  simula- 
tion is  typical  of  a medium-sized  advanced  fighter. 


plane 

(c)  No  pilot  v;ill  be  in  the  control  loop,  so  auto- 
matic operation  is  assumed 

(d)  The  variations  in  speed  of  sound,  mass  of  the 
aircraft,  and  atmospheric  density  are  negligible 

(e)  The  center-of-gravity  of  the  aircraft  does  not 
change  during  flight 

(f)  Perfect  navigation  data  on  the  state  of  the 
aircraft  is  used  (Ref  9:102) 
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II • Description  of  the  Control  System 


To  adequately  begin  discussion  of  the  proposed  ter- 
rain follov/ing  controller,  it  is  necessary  to  describe 
the  basic  components  of  the  entire  system.  The  remainder 
of  this  report  v/ill  further  explain  the  mechanization  of 
each  of  the  components.  Figures  la,  lb  and  Ic  illustrate 
the  control  system. 

First  of  all,  it  will  be  assumed  that  the  aircraft 
is  equipped  with  some  modern  sensor  device  which  provides 
accurate  range  and  altitude  data  for  the  terrain.  This 
can  be  achieved  by  using  a forward  looking  radar,  a laser 
ranger  or  even  a map-matching  terrain  data  processor. 

The  method  of  data  acquisition  is  not  important  to  this 
report,  as  long  as  good  terrain  data  is  available. 

Next,  this  reliable  range  and  altitude  data  is  fed 
into  the  reference  path  generator.  It  is  here  that  the 
reference  flight  path  is  constructed  in  the  range  domain. 
The  range  domain  is  chosen  in  order  to  more  directly  fit 
and  control  the  flight  path  with  respect  to  the  terrain. 

It  is  important  that  the  reference  path  is  constructed  to 
inherently  satisfy  any  acceleration,  flight-path  angle, 
and  clearance  constraints  imposed  upon  the  vehicle. 

Linear  variations  about  this  reference  trajectory  will  be 
considered  later,  and  any  limits  encountered  by  these 
perturbations  will  destroy  the  linear  behavior  and  degrade 
the  controller  performance.  The  reference  path  descrip- 
tion consists  of  the  path  altitude,  h,  and  its  first  four 
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(spatial)  derivatives  with  respect  to  range,  h , h , h 

IV 

and  h 

These  spatial  variables  are  used  for  all  subsequent 
reference  model  computations  to  obtain  the  desired  veloc- 
ity, the  reference  control  and  the  reference  state  vari- 
ables. The  velocity  generator  computes  the  desired 
velocity,  V^,  along  the  reference  path  by  using  some  form 
of  energy  management  criteria.  This  desired  velocity  can 
account  for  the  operational  speed  limits  of  the  aircraft 
as  well  as  the  desired  average  groundspeed  to  successfully 
meet  the  check  point  times  on  the  lov;  level  mission.  The 
desired  velocity  and  all  the  calculations  of  flight  path 
are  used  to  compute  the  complete  set  of  reference  states, 
x„.  These  reference  states  are  then  used  to  compute  the 
estimates  of  the  reference  controls,  Uj^,  needed  to  produce 
such  reference  states.  This  completes  the  generation  of 
all  the  necessary  reference  data  to  begin  the  simulation 
sequence. 

The  nonlinear  aircraft  differential  equations  are 
processed  to  simulate  the  actual  state  representation  for 
the  aircraft.  These  states  are  denoted  as  In  reality 

the  true  aircraft  states  cannot  be  determined  exactly. 

The  actual  states  would  be  estimated  from  a number  of 
measurements.  This  v;ould  involve  filtering  techniques 
and  various  sensing  devices.  However,  in  this  thesis 
the  control  system  errors,  rather  than  the  sensor  errors, 
are  of  major  interest;  therefore,  it  is  assumed  that 
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integration  of  the  differential  equations  of  motion  pro- 
vides the  exact  states  of  aircraft  motion  and  estimation 
theory  is  excluded. 

A comparison  is  next  made  between  the  actual  aircraft 
state,  x^,  and  the  corresponding  reference  states,  Xj^,  to 
compute  the  state  deviation  ^x.  These  state  deviations 
are  kept  small  by  the  selection  of  the  proper  feedback 
gain  matrix,  G,  More  explanation  will  be  given  to  this 
topic  in  Chapter  III. 

The  feedback  gain  matrix,  G,  is  calculated  through 
optimal  control  laws  (explained  in  Chapter  III)  and  the 
Ricatti  equation.  The  matrix,  G,  can  be  maintained  as  a 
constant  for  most  of  the  terrain  following  missions  to 
be  flown.  This  is  consistent  with  the  system  design  since 
the  state  deviations  are  small,  and  variables  such  as 
atmospheric  density  and  aircraft  mass  are  assumed  constant 
for  the  terrain  following  profiles  to  be  used  in  this 
thesis  (Gain  scheduling  is  possible  for  extended  missions). 

Once  the  feedback  gain  matrix,  G,  is  calculated  a 
feedback  control  signal,  ^u,  can  be  computed.  This  signal, 
^u,  is  achieved  by  multiplying  the  feedback  gain  matrix, 

G,  with  the  state  deviation,  ^x,  in  the  following  manner 
(Ref  4:151  and  13:194) 

^u  = G ^x  (1) 

The  feedback  control  signal,  ^u,  can  then  be  added  to  the 
previously  computed  reference  control,  Uj^,  to  provide  the 
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actual  control,  u.  This  control,  u,  provides  the  neces- 
sary commands  to  the  airv.raft  engine  and  aerodynamic 
control  system  to  maintain  the  aircraft  close  to  the 
reference  trajectory. 

With  the  computation  of  the  actual  control,  u,  the 
system  loop  in  Fig.  Ic  is  essentially  completed,  Nev/  data 
is  sensed  ahead  of  the  aircraft,  new  reference  parameters 
are  computed,  and  new  controls  are  activated  in  a sequence 
o.f  data  frames.  This,  basically,  is  the  overall  operation 
of  the  proposed  terrain  following  controller.  The  remain- 
der of  the  thesis  describes  the  components  of  this  control- 
ler and  their  interactions  in  more  detail.  The  results 
of  the  simulation  are  also  discussed. 
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III.  Optimization  of  the  Linear  Controller 


The  entire  approach  to  the  optimal  terrain  following 
controller  is  based  on  the  methods  of  modern  control 
theory.  This  theory  forms  the  foundation  for  the  prelimi- 
nary study  of  the  proposed  controller. 

The  Regulator  Problem 

The  function  of  a state  regulator  is  to  maintain  the 
states  of  a system  within  an  acceptable  deviation  from 
the  zero  reference  condition  by  using  required  control. 
This  is  the  type  of  feedback  controller  to  be  designed 
and  used  in  this  project. 

The  system,  in  this  case,  consists  of  the  deviations 
from  a time-variant  reference  model  and  thus  is  a linear- 
ized system, 

=■  A + B ^u.(+)  (2) 

In  Eq  (2) 

^x(t)  is  the  nxl  vector  of  perturbed  states 

^u(t)  is  the  mxl  vector  of  perturbed  controls 

A is  the  nxn  dynamics  matrix 

B is  the  nxm  control  distribution  matrix 

Both  A and  B are  time-variant,  but  they  are  assumed  to 
vary  slowly  to  allow  quasi-steady  state  solution  of  the 
controller  gains. 

The  reference  conditions  are  generated  by  the  path 
generator,  the  state  generator  and  the  control  generator 
as  described  in  Chapter  II.  In  addition,  the  reference 
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velocity  is  calculated  continuoualy  along  the  flight  path. 
For  this  thesis  the  initial  reference  velocity  is  647,3 
ft/sec  v/hich  is  Mach  .53  at  sea  level  conditions.  This 
initial  reference  velocity  can  be  used  as  a norninal  flight 
condition.  Since  variations  in  flight  condition  are  not 
large,  the  aerodynamic  stability  derivatives  are  assumed 
constant  in  the  region  of  interest  while  the  reference 
states  change  v/ith  time. 

The  perturbation  state  vector.,  , can  be  written  as 


where 

is  the  actual  state  vector 

and 

is  the  reference  state  vector 

and  the  perturbation  control  vector,  i'u,  can  be  written 
similarly  as 


where 

u is  the  actual  control  vector 

and 

is  the  reference  control  vector 

The  two  quantities,  Sx  and  represent  the  deviations 

from  their  respective  reference  conditions.  And  these  arc 
the  values  which  will  be  regulated  by  the  optimal  feedback 
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controller  developed  tor  this  thesis. 

For  ttie  Oi't  iinal  terrain  follov/inc  rennlator  an 
appropriate  cost  function  must  bo  dosipned.  The  most 
convenient  and  workable  form  is  tlie  quadratic  which 
weinhts  both  the  state  and  the  control  deviations  from  a 
reference  condition.  The  conventional  quadratic  cost 
function  v.'ill  be  (iNOf  '1:119) 

c 

This  form  sho'os  that  the  state  doA'iat. ions  and  the  control 
deviations  arc  directly  regulated  dependitig  on  the  form 
of  the  weighting  matrices  P and  R.  But  for  this  control- 
Icr,  not  all  of  the  state  deviations  are  to  be  regulated 
directly.  Instead,  an  o\.itput  deviation,  , will  be 
regulated  directly.  As  will  bo  seen  later,  this  merely 
represents  a particular  choice  of  the  P matrix. 

The  output  is  defined  as 

in  C (6) 

/ ^ 


where  C is  the  transformation  matrix.  Using  the  output 
deviation,  Sy , the  cost  function  can  be  written  as 


J 


R cit 


o 

The  selection  of  the  two  weighting  matrices,  Q and  R 
be  described  at  the  end  of  this  chapter. 

The  transpose  of  the  output  vector,  ^y,  is 


(7) 

, will 
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(B) 


S/-  (cu)  = 5/c^ 

Substitution  of  Sq  (G)  and  Gq  (6)  into  Eq  ( "y ) produces  the 

revised  cost  function 

T-t>^ 

J - “ QC  Sx  +W^R$U^c!i‘  (9) 

o 

Comparison  of  Eqs  ([>)  and  (9)  shov/s  that  they  arc  oquiva- 
T 

lent  if  P = C QC.  The  cost  function  in  Eq  (9)  is  also 
in  the  form  to  use  the  Ricatti  equation  to  determine  the 
Ricatti  gain  matrix,  S. 

The^  Ricatti  Equation 

The  solution  to  the  regulator  problem  requires  solu- 
tion of  the  Ricatti  equation  (Ref  4:167).  For  simplicity, 
it  is  assumed  that  the  dynamics  matrix,  A,  and  the  control 
distribution  matrix,  B,  are  constant  matrices.  This 
assumption  is  not  strictly  true  and  results  in  a sub- 
optirnal  sot  of  gains.  In  the  cost  function  of  Eq  (9)  the 

weighting  matrix,  R,  is  assumOu  constant  as  is  the  product 
T 

matrix,  C QC.  'With  those  assumptions  the  stead,,  state 
Ricatti  equation  is  written  in  the  form 

S = 0 = -SA-a’s  + 9BR'B^S-c"0C  do) 

This  equation  above  implies  that  the  Ricatti  gain  matrix, 
S,  is  constant  for  the  optimal  terrain  I'ollowing  regulator 
(Ref  4:167),  If  terrain  follovjing  occurred  for  prolonged 
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periods,  the  gains  could  be  scheduled  for  varying  flight: 


conditions.  The  solution  of  Eg  (10)  is  accomplished  by 
using  the  OPTSYS  program  described  in  Appendix  E, 


Linear  State  Feedback 

V/hen  solving  for  the  feedback  gain  matrix,  G,  from 
Chapter  II,  the  linear  control  lavj  of  Eq  (1)  can  be 
v;ritten  in  vector  form  as 


r:  G $ 


(11) 


Thrs  value  for  can  be  substituted  in  Eq  (2)  to  obtain 

~ A + BG-^x  112) 


Factoring  out  fruia  Eq  (12)  yields 

^ X — ^ A + B $ X (13) 

The  term,  A+BG,  :s  knov;n  as  the  closed  loop  dynamics 
matrix  because  this  term  incorporates  the  feedback  gain 
matrix,  G,  as  seen  in  Fig,  1, 

A connection  may  nov;  be  made  between  the  feedback 
gain  matrix,  G,  and  the  Ricatti  gain  matrix,  S.  Since  S 
is  the  steady-state  Ricatti  gain  matrix  it  can  be  shov;n 
that  (Ref  4:167  and  13:221) 

G'=-r'b's 


and  the  steady-state  linear  control  lav/  of  Eq  (il) 
becoruos 


Si  - -r’b^sSx 


(15) 


Thin  equation  above  forinn  the  banic  conti'ol  law  for  the 
linear  requlator  problem  considered  in  this  thesis.  It 
is  llq  (15)  which  updates  the  cotitrol  for  the  terrain 
following  regulator  according  to 
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no  t.ei'm  i nn  t i on  of  0 and  11  t 'in tici  cos 

Ti’iero  is  no  one  standard  method  for  deton, lining  the 
exact  components  of  tlic  Q and  R weighting  matrices  for 
any  given  jn-oblom.  It  requires  some  itc'rations  and 
educated  guesswork  to  determine  the  best  overall  wclgnts. 
Both  the  Q and  R matrices  v/jll  b(.'  used  in  diagonal  form 
only,  'ihis  is  done  because  normally  only  tt'io  main  diago- 
nal Clements  arc  significant  (Ref  A:l‘\0  and  10:33).  The 
Q matrix  weights  cacti  dcviati'>n  in  tlic  four  output  vari- 
ables and  tne  R matrix  v;oights  the  deviation  in  the  two 
control  variables. 

The  method  of  Bryson  and  llo  (Ref  d:ld9)  can  be  used 
\o  obtain  vjorking  values  for  ilie  diagonal  elements  of  both 
Q and  U.  This  mottiod  utllir.es  the  setting  of  maximum 
accoiitablc  values  on  each  cleinont  of  tlio  output  deviation 
vector,  £ V , and  tlic  control  deviation  vector,  £_u.  Corre- 
sponding elements  of  the  v/ciytiting  matrices  can  tlicn  be 
used  in  the  cor.t  function,  0.  Tlio  equation.':  for  dotcr- 


18 


minincj  Lhe  diayonaJ  element  weiglits  are 
1/Q^j  a the  nuiximuin  acceptable  value  of 


and 

l/Rf^  the  maximum  acceptable  value  of 


( 10) 


It  should  be  noted  that  these  tv;o  equations  will  produce 
a ’'first  guess"  set  of  diagonal  elements  v^ihich  can  be 
adiu^ted  to  obtain  rcasona'ulo  performance  for  the  cost 
function.  Once  the::c  wcigliting  matrices  liave  been  esta- 
blished tlie  resulting  Hicatti  gain  matrix,  S,  is  auto- 
matically fixed  accot'diny  to  the  Kicatti  equation,  llq 
(10).  This  matrix,  S,  in  turn,  directly  affects  the 
overall  value  of  the  feedback  gain  matrix,  0,  as  seen  in 
Eq  (M).  Therefore,  any  change  in  the  diagonal  elements 
of  either  0 or  K will  create  multiple  changes  in  the 
subsequent  gain  matrices  of  the  system. 

The  reason  bryson  and  Ho  oquation^•.  co,n  bo  corisidered 
a "first  guess"  estimate  is  that  certain  selected  Q and  R 
diagonal  elements  may  result  iri  control  values  that  arc 
not  realisable.  In  this  fashion  the  Q and  K matrices  must 
be  judiciously  adjusted  to  obtain  good  perfor'rnance  without 
adversely  aliecting  tt>c  flying  qualities  or  the  realisable 
control  deviations. 

If  tlio  diagonal  elements  of  Q are  increased  the 
controller  tetujs  to  become  strongi'r  and  restrict  the 
system  deviations  closer  to  ccro.  Tl\is  can  also  be 
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achieved  by  dccrcaaing  the  element  values  of  R , A 
stronger  terrain  following  regulator  will  result  in  a 
decrease  in  the  response  time  of  the  overall  system.  But 
this  is  usually  accomplished  v;ith  a corresponding  increase 
in  peak  overshoot.  So  these  tradeoffs  must  be  considered 
when  designing  the  best  values  for  the  diagonal  elements 
of  the  Q and  H weighting  matrices. 
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IV 


Tho  f>tat:e  Vnrinhlos  and  Kquatlonr,  of  MoV.  ion 


Tho  purpose  of  this  chapter  will  be  to  denote  the 
variables  comprising  tho  state  vector  and  tht  control 
vector  used  in  this  thesis.  Additionally,  this  chapter 
will  describe  the  linearised  equations  of  motion,  A full 
development  of  tho  equations  of  motion  will  not  be  pre- 
sented here,  but  the  reader  is  directed  to  References  6, 
7,  12,  14  and  16  if  a fuller  explanation  is  needed. 

Only  the  longitudinal  equations  of  motion  will  be 
used  in  this  controller  design  since  tho  longitudinal 
dynamics  are  the  most  important  in  a terrain  following 
flight  profile.  The  flight  regime  V'/ill  be  restricted  to 
the  x-z  plane  of  motion  as  seen  in  Fig,  2.  Lateral  con- 
trol is  usually  termed  terrain  avoidance  and  is  beyond 
the  scope  of  this  thesis.  Omission  of  the  lateral  equa- 
tions of  motion  will  not  degrade  the  controller  design 
nor  provide  inaccurate  results.  And  because  of  near- 


Figure  2.  Two  Dimensional  Terrain  Following  Profile 
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symmetry  of  the  aircraft  the  longitudinal  equations  of 
motion  may  be  separated  from  the  lateral  equations  of 
motion  (Ref  17:218). 


Assumptions  and  Coordinate  Axis 

While  developing  the  longitudinal  equations  of  motion 
the  aircraft  is  assumed  to  be  a rigid  body.  The  variation 
of  the  pitch  angle,  , is  arsumed  small  so  that 

COS  S©  = 1 

and 


i 

I 


t 


I .1 


1 

i 
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''tcViiib  jin.-M.  n-  ^ | 


And  the  forward  velocity,  V,  is  approximated  by 


M - V\F+  ^ V 

The  coordinate  axes  to  be  used  will  be  a set  of 
stability  axes  that  is  fixed  with  respect  to  the  reference 
trajectory  velocity  vector.  The  relation  to  the  x-z 
body  axes  is  shown  in  Fig,  3. 

State  Variables 

For  the  design  of  this  terrain  following  controller 
there  are  seven  state  variables  chosen  to  represent  the 
states  of  the  system.  These  are 

U the  longitudinal  forward  velocity 

w the  normal  velocity 

q the  pitch  rate 

0 the  pitch  angle 

£ the  elevator  deflection  angle 

A the  thrust 

h the  altitude 

These  seven  state  variables  will  be  used  throughout  this 
thesis  and  will  be  collectively  designated  as  the  vector 

21* 

Two  of  the  state  variables  will  be  directly  control- 
led in  this  design  project.  These  two  command  values 
are  the  controls 
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the  commanded  clovutor  deflection  angle 
the  commanded  thrust 


and  can  be  denoted  as  the  control  vector  u. 


Perturbation  Theory 

For  the  purpose  of  designing  the  tracking  controller, 
the  nonlinear  longitudinal  equations  can  be  linearized 
about  a particular  reference  trajectory  and  the  motion 
described  in  terms  of  the  variations  in  seven  state 
variables.  This  is  illustrated  in  Fig.  4. 


Figure  4.  Perturbation  about  a Reference  Trajectory 
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The  use  of  small-disturbance  theory  car,  be  used  to  pro- 
duce the  linearized  perturbation  equations  of  motion 
(Ref  7; 120).  All  the  perturbation  quantities  and  their 
respective  derivatives  are  assumed  small  so  their  squares 
and  products  are  negligible  when  compared  to  the  first- 
order  linearisation  terms.  A linearised  perturbation  in 
angle  of  attack,  , is  illustrated  in  Fig.  5. 


Figure  5.  Linear  Perturbation  in  Anglo  of 
Attack 


Linearised  Perturbation  Equations 

The  linearized  perturbation  equations  for  each  of 
the  seven  state  variables  can  now  be  written  using  dimen- 
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siona]  stability  dorivatlvos.  These  dimensional  stability 
derivatives  are  explained  in  Appendix  B. 

The  fact  that  the  equations  are  linearized  about  an 
arbitrary  trajectory,  rather  than  straight  and  level 
flight  as  is  frequently  done,  means  that  the  coefficients 
in  the  linearized  equations  would  depend  on  the  rofcrenco 
trajf  story  states.  This  implies  that  the  coefficients 
would  vary  with  time.  In  order  to  develop  a controller 
with  constant  gains  the  coefficients  must  be  constant. 

The  assumption  required  for  this  is  that  the  reference 
flight  path  angle,  » ^nd  angle  of  attack,  ^irc 

small  enough  that  their  producLs  v;ith  a perturbation  state 
can  be  considered  of  second  order.  In  this  v;ay  the 
perturbation  equations  raduco  to  tlie  usual  form.  Therefore, 
for  the  purpose  of  this  study  the  usual  linearization 
about  straight  and  level  flight  was  considered  accurate 
enough  to  use  for  the  linear  perturbation  model. 

For  tlie  perturbation  in  longitudinal  velocity,  U: 

Sv  = 't-  Xv^^VsJ  + ^ ( 19) 

For  the  perturbation  in  normal  velocity,  w: 

v;here  U is  the  nominal  forv.ard  velocity.  For  the 
o 

perturbation  in  pitch  rate,  q: 
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(22) 


For  the  perturbation  in  pitch  angle,  © : 

Se  = % 

For  the  perturbation  in  elevator  deflection  angle, 


JS  = - SS,)  ( 

where  'T  is  the  time  constant  associated  with  the 

* c 

elevator  control,  arid  is  the  elevator  control 

command  based  on  the  feedback  of  perturbation  states. 
For  the  perturbation  in  thrust,  ^ : 


(2'i) 


where  *T  is  the  time  constant  associated  witJi  tfie  thrust 
P 

control,  A » 3nd  is  the  perturbation  t’-.rust  command. 

For  the  perturbation  in  altitude,  h: 

12  -t  (25) 

These  equations  above,  Eq  (19)  ~ Eq  (25),  constitute 
the  long itudir.al  perturbation  equations  of  motion  in 
dimensional  derivative  form.  Those  equations  can  be 
rewritten  in  the  matrix  form  of  the  linear  differential 
equation  as 

= A Sx  -V-  BW  (26) 

where  A is  the  open  loop  dynamics  matrix.  Gathering 
Eq  (19)  through  Eq  (25)  into  matrix  form  results  in: 
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A= 


X,  -g  Xj  X,  0 


M,.,  Mo 


0 0 


0-10 


0 Z;  0 

0 Mj  0 0 

0 0 0 0 
O -Zr,  0 O 

O 0 'Yr,  0 

U 0 0 0 


(27) 


Forcinq  Function 


The  two  aircraft  control  commands,  and  , 

act  as  forcing  functions  for  the  linearised  system  of 
equations.  This  forcing  function  consists  of  a control 
distribution  matrix,  B,  multiplied  by  the  control  command 


vector.  The  matrix  U takes  the  form  of 
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(28) 


The  commanded  controls  must  also  be  v/ritton  in  the  pertur- 
bation form  as 


(29) 


This  completely  describes  the  terms  in  perturbation 
equation  (26)  which  is  identical  to  Eq  (2)  in  Chanter  III, 
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^ * Cubic  :Iplino  Ivoferoncn  Patli 

This  chapter  discusser,  the  main  points  x’egardincj  the 
generation  of  the  optimised  reference  flight  path.  It  is 
this  reference  path  v.'hich  the  optimal  terrain  following 
controller  is  designed  to  follow  as  closely  as  possible. 
All  of  the  information  in  this  chapter  originates  in 
Referei'ice  7 since  the  reference  path  generation  is  the 
essence  of  Major  Funk’s  doctoral  dissertation.  As  a 
result,  this  chapter  v.'ill  not  detail  the  reference  path 
generation,  but  will  highlight  those  main  ideas  of 
importance  to  the  aesign  of  an  optimal  terrain  following 
controller. 


Cubic.  Splines 

A cubic  spline  function  consists  of  a sequence  of 
cubic,  polynomial  segments  v/hich  together  form  a continuous 
function  having  continuous  first  and  second  derivatives 
(Ref  liS-l).  Mathematical  spline  functions  are  ideally 
suited  to  digital  trajectory  computation.  They  provide 
smooth  continuous  curves  that  can  bo  specified  uniquely 
by  a discrete  sot  of  parameters. 


A cubic  polynomial  can  bo  denoted  as 


for  i ^ 1,  2,  n*"! 


where 


X 
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A typical  cubic  cplinc  function  is  shown  in  Ficj,  6, 


f i<;uro  6.  Cubic  Spline  Function 


The  function  points  y^)  called  knots. 

The  first  and  second  derivatives  can  be  obtained  by 
dif f orenti-ating  Eq  (30).  For  the  interior  knots  (x^ 
through  x^  in  Fig  , G)  the  first  and  second  derwatives 
must  be-  identical  for  each  adjoining  interval-.  This 
means  the  derivatives  at  knot  <X2,  y2)  must  bo  the  same 
for  to  X2  ana  for  interval  x-,  to  x^.  This  insures 
continuity.  Therefore,  the  three  unknovfns  for  each 
segment,  A^ , and  C^,  can  be  computed  by  enforcing  the 

contiriuity  conditions  at  each  interior  point  and  by 
specifying  the  slope  at  each  endpoint.  In  this  manner  a 
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ur.iquo  cubi...  spline  function  can  be  constructed  to  ropro- 
sent  the  (jlven  data  as  a smooth,  continuous  function. 

For  a terrain  following  profile  it  is  important  to 
let  the  knots  coincide  with  the  terrain  data.  This  proce- 
dure insures  that  the  cubic  function  passes  through  each 
of  the  data  points. 

Tc." r rn  n Data  and  Clearance  Curve 

The  sampled  terrain  data  is  defined  in  the  range 
domain  rather  than  the  time  domain.  This  is  ideal  for  a 
terraifi  follov;ing  flight  profile  since  terrain  elevation 
is  naturally  described  as  a function  of  range  rather  than 
time.  Therefore,  Fig.  6 can  represent  the  known  terrain 
data  if  the  ordinate  becomes  terrain  elevation  and  the 
abscissa  is  range.  This  is  illustrated  in  Fig,  7, 


Figure  7.  Sampled  Terrain  Data 
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The  reference  path  generation  bogina  ^ith  the  input 
of  sampled  terrain  data.  These  data  point-  consist  of 
th(5  elevation  of  the  terrain  at  specified  range  intervals. 
For  instance,  the  elevation  (in  feet)  is  specified  for 
each  2000  feet  of  range  along  the  flight  profile.  The 
reference  path  generator  algorithm  normalizes  the  raw 
data  and  constructs  a cubic  spline  function  which  passes 
through  the  given  set  of  terrain  points.  This  is  called 
the  terrain  curve.  Next,  the  minimum  amount  of  clearance 
above  the  terrain  curve  is  specified.  This  designated 
clearance  is  used  for  determination  of  the  cubic  spline 
c].carance  curve.  Both  the  terrain  curve  and  the  clearance 
curve  arc  shown  in  Fig,  8 below. 


Figure  8.  Terrain  and  Clearance  Curves 
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The  altitude  along  the  cubic  spline  clearance  curve 
is  designated  h,  which  is  a function  of  range.  Since  the 
clearance  curve  is  processed  as  a cubic  spline  path,  h 
v/ill  represent  a cubic  spline  function.  The  slope  of  tlie 
clearance  curve,  h,  is  defined  as  the  first  derivative  of 
h with  respect  to  range  (spatial  derivative)  and  desig- 
nated as  either  h'  or  s. 


w = 


iL 

JR 


- s 


The  second  spatial  derivative  if 


JR 


I = k 


whore  h or  k is  called  the  curvature. 
And  the  third  spatial  derivative  is 


r _ A 


- P 


lU 


(31) 


(32) 


(33) 


where  p or  h is  called  the  kink.  As  mentioned  in  the 
previous  section  on  cubic  spline  functions,  the  variables 
h,  s and  k are  continuous  and  p is  piccev/ise  continuous. 


The  Reference  Path  Optimal  Control  Prob 1 om 

The  cubic  spline  clearance  curve  can  be  considered 
a minimum  clearance  curve  over  the  terrain.  It  would  be 
extremely  difficult,  if  not  impossible,  for  a manned  air- 
craft to  track  this  minimum  clearance  curve  efficiently 
within  the  aircraft  performance  limitations.  Therefore, 
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the  cubic  spline  clearance  curve  must  be  optimized  by 
minimizing  a performance  index  using  optimal  control 
theory.  Por  this  optimal  control  problem  it  is  possible 
to  specify  certain  constraints  on  k,  s and  p so  that  the 
limits  of  the  aircraft  are  not  exceeded  and  the  aircraft 
does  not  descend  below  the  minimum  clearance  curve.  The 
performance  index  rainiinizes  the  deviation  from  the  clear- 
ance curve  while  enforcing  the  specified  constraints  along 
the  curve. 

T'ne  reference  path  optimization  problem  is  solved 
subject  to  the  follov/ing  inequality  constraints  along  the 
entire  range  (Ref  9:35) 


s , « s ^ 

min  “■ 


^min  - ^ - 


P • ^ p ^ 

*^min  — ^ 


“^max 

max 

D 

Mnax 


The  limits  on  k are  determined  from  the  normal  accelera- 
tion limits  and  the  fiominal  aircraft  speed  by  the 
following  equation 


k 


This  equation 
small.  The  k 


is  accurate 
limits  tend 


(34) 


as  long  as  the  slope,  s,  is 
to  occur  at  the  tops  of  peaks 


and  near  the  bottoms  of  valleys  v/here  the  slope  is  nearly 
zero.  So  Eq  (34)  is  a good  approximation. 


/ 

Framing  of  the  Optimal  Control  Problem 

It  would  be  very  impractical  to  solve  an  extremely 
large  optimal  control  problem  for  the  entire  terrain  fol- 
lov/ing  flight  profile.  As  a result,  the  optimal  control 
problem  must  be  broken  into  specified  segments  called 
frames  (Ref  9:44).  ' 

The  framing  procedure  is  illustrated  in  Fig.  9 v/here 
there  is  no  overlap  betv;een  frames.  For  Frame  1,  essen- 
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Frame  3 


I Frame  2 ^ 


Frame  1 


% 


K 


Figure  9.  Framing  v;ith  no  Overlap 


tially  an  optimal  control  problem  is  performed  over  the 
range  interva)  to  At  range  a new  optimal 

control  problem  is  begun  for  the  range  interval  for 
^1  ^2*  framing  procedure  with  no  overlap  v;ould 

v/ork  well  if  the  input  terrain  data  v;ore  perfect.  But 
because  of  radar  sliadowxng  the  present  state  of  the  art 


36 


. , • <>CTB-'W>WgWO>WftWWW^ 


does  not  allov;  the  accurate  solution  vijith  no  frame  over~ 
lap.  There  must  be  an  overlap  between  frames  so  that  as 
updated  terrain  data  ahead  of  the  aircraft  is  sensed,  a 
new  optimization  problem  can  be  solved  to  update  the 
previously  computed  flight  path  signals.  The  computer 
algorithm  developed  by  Major  Funk  optimizes  the  entire 
reference  flight  path  by  overlapping  frames  as  can  be 
seen  in  Fig.  10,  This  overlap  allov/s  for  a continuous 
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Figure  10.  One-half  Frame  Overlap 


update  of  the  reference  flight  path  as  new  terrain 
information  is  detected. 

Within  each  frame  a certain  number  of  control  points 
are  specified.  These  control  points  correspond  ho  the 
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points  at  which  tho  constraints  on  h,  s,  k and  p must  bo 
satisfied.  As  the  number  of  control  points  per  frame 
increases,  the  computation  time  also  increases.  So  there 
is  a tradeoff  betv;een  using  a large  number  of  control 
points  which  forces  the  reference  path  to  more  closely 
follow  the  terrain,  and  using  a small  number  of  control 
points  which  reduces  the  computation  time  and  the  control 
rate. 

In  addition  to  the  control  points,  there  are  numer- 
ous other  parameters  which  must  be  specified  in  each 
frame  of  the  solution  to  the  optimisation  problem.  These 
various  parameters  will  not  be  described  hero  since  they 
are  listed  in  Reference  9,  page  54.  3ut  the  availability 
of  such  a variety  of  parameters  permits  the  flexibility 
in  the  generation  of  the  optimized  cubic  spline  reference 
path  for  various  aircraft  and  missions.  This  reference 
path  generator  can  account  for  almost  all  the  practical 
constraints  for  any  type  of  flight  vehicle  on  varying 
teirain  following  profiles^ 

Output  of  the  Reference  Path  Generator 

The  output  of  the  reference  path  generator  consists 
of  the  optimized  cubic  spline  reference  path  as  shovm  in 
Fig.  11,  It  can  be  seen  that  the  reference  path  docs  not 
coincide  v;ith  the  minimum  cloarar.iie  curve.  The  main 
reason  for  the  difference  bet-;en  the  two  curves  is  that 
the  reference  path  must  satisfy  all  the  specified  con- 
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straints  along  the  flight  profile.  Generally  some  por- 
tions of  the  clearance  curve  v;ill  have  excessive  curva- 
tures or  slopes  which  would  exceed  the  acceleration  and 
flight-path  angle  constraints.  Whenever  the  clearance 
curve  satisfies  the  specified  constraints  the  two  curves 
will  coincide. 

The  actual  numerical  output  of  the  reference  path 
generator  consists  of  the  optimised  path  altitude,  h,  and 
its  sloi^e,  h'  , at  each  discrete  increment  of  range,  R^. 
This  path  data  is  then  input  to  an  interpolation  routine 
which  reconstructs  the  full  cubic  spline  using  a direct 
slope  determination  technique  (Ref  9;61  and  18:54).  This 
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cubic  spline  interpolation  routine  calculates  the  follow- 
ing parameters  at  any  desired  range  along  the  flight 
profile : 

(a)  h - altitude 

(b)  h*  - slope 

Cc)  h - curvature 

(d)  h**'  - kink 

These  four  quantities  can  then  be  used  in  the  nonlinear 
equations  which  represent  the  cubic  spline  terrain 
following  controller  outlined  in  Chapter  VI. 


Some  Important  Relationships 

Before  describing  the  nonlinear  reference  model  it 
will  be  helpful  tw  show  some  important  relationships 
between  the  time  and  spatial  derivatives  along  the 
reference  path.  Since  altitude,  h,  is  really  a function 
of  range  it  can  be  written 

V\  ~ (35) 

\>;here  range,  R,  is  a function  of  time.  From  Fig.  12  the 
time  derivatives  of  range  and  altitude  are  shov;n  in  rela- 
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tion  to  the  velocity,  V.  These  relationships  can  be 
written  as  (Ref  11:306  and  15:64) 


U = V 


and 

R =:  V cos  ^ 

But  h can  also  be  written  as 

K - iL  iE.  ::  K'V  cos'!! 

JR  Jt 


(36) 


(37) 


(38) 


which  is  another  representation  of  the  same  quantity  in 
Eq  ( 36 ) . 

/ 

The  time  derivative  of  slope,  h , can  bn  computed 


from 


d(t')  , Jk'  JR 
J+  ■ JR  i\ 


(,''Vcosl5 


(39) 


These  equations  above  will  be  important  for  the  nonlinear 
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This  chapter  describes  the  nonlinear  equations  used 
to  construct  the  reference  model  for  the  optimal  terrain 
following  controller.  In  Chapter  IV  the  linear  model  with 
the  perturbation  equations  was  presented.  But  the  non- 
linear equations  must  be  used  in  the  design  simulation 
since  the  nonlinear  model  more  accurately  represents  the 
actual  aircraft  flight  regime.  The  nonlinear  response 
and  stability  characteristics  provide  a more  realistic 
controller  operation.  The  purpose  of  this  chapter  is  to 
show  the  development  of  the  equations  leading  to  the 
reference  model  state  vector,  Xft»  truth  model, 

X.,  Tnen  these  two  state  vectors  can  be  used  to  obtain 

—A 

the  commanded  control  vector,  u,  necessary  for  tracking 
the  cubic  spline  reference  path.  A partial  block  diagram 
of  this  development  is  shown  in  Fig.  13  below. 


Figure  13.  Computation  of  Commanded  Control,  _u 
An  attempt  has  boon  made  to  organize  this  chapter  sequen- 
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tially  through  the  equation  derivations  in  the  order  in 
which  they  are  required. 

Reference  Path  Generation 

The  cubic  spline  reference  path  generator  described 
in  Chapter  V provides  the  necessary  altitude  and  deriva- 
tive data  to  construct  the  optimised  flight  path.  The 
path  generator  calculates  the  following  reference  quan- 
tities in  the  range  domain: 

(a)  - reference  altitude 

(b)  h^  - slope 

i/ 

(c)  h^  - curvature 

(d)  ” kink 

where  the  subscript,  d,  denotes  desired  or  reference 
quantity.  These  four  quantities  are  essential  to  the 
compvitation  of  many  other  reference  conditions. 

Reference  Velocity  Generator 

To  correlate  time  and  spatial  relations  along  the 
flight  path,  velocity  must  be  known.  The  desired  velocity 
along  the  reference  path  could  be  computed  from  energy 
management  considerations  once  the  path  is  known.  For 
simplicity,  a constant  energy  assumption  is  made  here. 

If  the  terrain  is  not  too  severe  and  the  vehicle  kinetic 
energy  is  large,  this  is  a reasonable  assvimption.  This 
allows  an  interchange  of  kinetic  and  potential  energy 
with  a corresponding  reduction  in  throttle  activity  over 
that  required  for  constant  speed  flight. 
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StaL'tir.q  with  tho  initial  velocity,  V , and  the 

-’o’ 

initial  altitude,  h^,  the  total  specific  energy  of  the 
system  is  coinputed  by  sumrninn  the  specific  potential  and 
kinetic  energies  as 

= Z 3^0  (40) 

The  cubic  spline  reference  path  generator  provides  a 
value  for  desired  altitude,  , at  any  point  along  range. 
So  the  desired  velocity,  , may  bo  computed  by  rearrang- 
ing the  constant  energy  equation 

V,  = V2(e, 

This  is  the  equation  for  calculating  desired  oi  reference 
velocity  along  range. 


Reference  Flioht  Path  Annie  and  First  Time  Derivative 

The  flight  path  angle,  , can  be  computed  using  the 
trigonometric  relation  shov;n  in  F'ig,  14. 
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Since  in  computed  by  the  rciorenco  poth  goncralo::  ttie 
do:;  1 red  t light,  path  angle,  ^ ^ , ir-  ' 
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bo  written 
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(•lO 


Wl  u;r«,: 
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R = V cos 


f 

An  expro.nr.ion  tor  ^ can  bo  derived  by  t.iKiny  ttie  r.pcitial 
dcrlvotivo  oi  t;otli  ride.';:  ol  Uq  (•'12)  » Thin  re.';wit:;  in 


If'  = K"  «$'"!( 


(•Ih) 


Sub:;litvii-.i  .Kj  hq  (.37)  and  hq  (-lb)  into  hq  (■•’•!)  yicldr  the 
ex-pre.'-.;;  iui'i  for  the  1 ii  ;'.l  time  d.-rivatlve  of  ^ 

\ •*  Xi  cos"^ 


(•lb) 


Time  Per  ivah  ivc:'.  oi'  Voioc-i<:v 

Now  that  the  flight  path  angle  and  its  first  time 
derivative  arc  derived,  it  is  possible  to  obtain  the  first 
and  sccoiKi  time  derivatives  of  desired  velocity. 

The  assumption  of  constant  energy  implies 


E = O 

o 


(47) 


After  differentiation  of  both  sides  of  Eq  (40)  the  first 
time  derivative  of  velocity  becomes 


<1,  = 


(48) 


From  Eq  (38)  in  Chapter  V the  h,  term  may  be  replaced  by 


(49) 


Substituting  Eq  (49)  into  Eq  (40)  yield: 


\ -- 


3 '■'i 


(50) 


This  is  the  expression  for  the  first  time  derivative  of 
desired  velocity.  The  second  time  derivative  can  bo 
derived  in  a similar  manner. 


If  Eq  {^2)  is  written  with  the  reference  sul:)script 


d,  it  bocornos 


hj  = +<v«  \ 


(51) 
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Thi::.  form  of  pcrmitG  r.ubstitution  in  Eq  (50)  to  pro- 
duce 


(52) 


Differentiating  both  sides  of  Eq  (52)  v/ith  respect  to 
time  results  in 


% = -9 

« 

Substitution  of  Eq  (46)  for  shows  that 

v,= 


(53) 


(54) 


Eq  (54)  is  the  equation  for  the  second  time  derivative 
of  desired  velocity  and  this  equation  v/ill  be  used  in 
the  nonlinear  model. 


Definition  of  K and  itr.  Dor ivat Ivor. 

Since  one  important  term  reappears  in  many  subse- 
quent equations  it  is  convenient  to  define  the  quantity, 

K,  called  the  mattiemat ical  path  curvature  (it  is  inversely 
proportional  to  the  instantaneous  radius  of  curvature 
(Ref  9:149)). 


i/  ^ I *'  ^ V 

K = cos 


(55) 


The  first  derivative  of  K with  respect  to  range  is 


/ 


K 


<//  ^ 7,  ^ 

Uj  cos 
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Substitution  of 


r 

Eq  (<15)  for  yields 

K ~ COS^  ” 3(^1  ) C05  )(j  S\ri 


But  Eq  (57)  can  bo  written  in  terms  of  K as 


K’  = cos^  'll  " Sm'Jj 


(57) 


(58) 


Eq  (50)  is  the  expression  for  the  spatial  derivative  of 
K,  (In  subsequent  equation  derivations  the  subscript,  d, 
vJill  be  omitted,  but  it  is  still  assumed  tliat  all  quan- 
tities will  be  those  desired  for  the  reference  model). 

The  time  derivatives  of  K arc  now  obtained  from 
differentiation  of  Eq  (55)  v^ith  respect  to  time  so  that 


K = 4l!i}  coi)i 

At 

But  it  can  be  seen  that 

jCC)  _ JR, 
J+  " JR  at 


- is 


(59) 


(60) 


v/hich  equal; 


a(u") 

at 


Hf 


“ U Vcos'^ 


(61) 


If  the  V ccuation,  Eq  (46),  is  written  in  terms  of  K 
it  becomes 


i r.  KV 


(62) 
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Substituting  Eq  (Gi)  and  Eq  (62)  into  Eq  (59)  shows  that 


the  first  titrie  derivative  of  K is 

tti 


K = (>"Vcos^)!  -5W'KVsm1S  cos'li 

And  this  can  be  written  more  compactly  as 

k = V(k'co5‘'x  - 3WK'’) 


(63) 


(64) 


The  second  time  derivative  of  K results  from  taking 
the  derivative  of  both  sides  of  Eq  (C4). 

k = V(k'co6^  -3WK')  +’((CVcos"^  - 
Sirv-iCi  -3k''Vco55!K'-4U'Kk') 


(G5) 


v;here  h =0  for  a cubic  spline  function. 

After  much  manipulation  and  rearrangement  of  Eq  (G5) 
it  can  be  rov/rittcn  in  slightly  shorter  form  as 

(GG) 

Second  Time  PGrivativo  of  Fliqnt  Path  Anglo 

v;ith  all  the  equations  for  K it  is  now  possible  to 
derive  the  second  time  derivative  of  the  fliglit  path 
angle.  Differentiation  of  Eq  (G2)  v/ith  respect  to  time 
is 


k =k(^-iovkV^^ 

u'w/x-skfu 

sk'l 

\V  / 

> • 

V 

/ 

H - KV  + V K 

Substitution  of  Eq  ( G4 ) into  Eq  (G7)  produces 


(67) 
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(63) 


^ ^ (u  COS  ^ + VK 

This  equation  is  the  expression  for  the  second  time 
derivative  of  and  is  used  later  for  computation  of  the 
angular  rates  for  uhc  nonlinear  equations. 

Forces  in  the  Lift  Direction 

This  section  derives  one  important  equation  v/hich 
v/ill  bo  used  for  ott'ior  reference  model  computations. 

This  equation  results  from  the  summation  of  all  forces 
in  the  lift  direction.  From  the  forces  shov/n  in  Pig.  15 


the  force  equation  is 


|YV  L "t  T S\r\  OC 


COS 


(G9) 
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c.  + 


sir^oc  = 3 


( 73) 


1-  t 

The  lift  coefficient,  C^,  con  be  expressed  as 


C “ C,  -V-  C,  <><■+•  C ^ + C A<1  (74) 


v;horo  fl  denotes  Mach  number.  Using  the  small  angle 
approximation  and  substituting  Eq  (74)  into  Eq  (73) 
yields 


C,  (75) 

Eq  (75)  is  the  first  of  two  equations  v/hich  v/ill  be  solved 
as  a simultaneous  sot  later  in  this  chapter. 


Estimated  Angle  of  Attack  and  Per ivativor.; 

In  order  to  obtain  an  estimate  of  the  angle  of 
attack  (call  this  o<j^)  LSio  following  conditions  are 


c o 


on  Eq  (75) 


c 


L 


o 


This  is  done  so  the  resulting  expression  for  the  estimated 
angle  of  attack,  c><^,  can  be  used  for  computing  the  time 
derivatives  of  o(  needed  in  subsequent  equations  (the 
actual  reference  value  for  o<.  is  obtained  later  v/ithout 
assuming  the  zero  conditions  above).  So,  after  the  scro 
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conditions  arc  sot  in  Cq  (75)  the  resulting  expression 


tor  estimated  anqle  of  attack  is 


^ oi  /VI 


(76) 


At  this  point  it  is  convenient  to  define  a variable, 


c^.  = = 


Ca  = 


+ — 
m ^ '■« 


-I 


(77) 


Using  Eq  (77)  t’no  expression  for  o<  becomes 


(78) 


Differentiation  of  Eq  (78)  with  respect  to  time  results  in 

«E=  Ca(k\;V5cos^')+^a('<'^^+2KVV'- 

^ sin')}  'i  j C79) 


An  expression  for  ^ 


"A 


Ca  = - 


which  can  be  written 


^SQ„  + 7 


tn 


(00) 


/VI 


= -Ca^VVSCu 


‘«C 


(81) 


Substituting  Eq  (81)  into  Eq  (79)  yields 


U =.  /KM^+gcosVj  + 


+2KVV  -qsv.\'«l! 


(82) 


But  it  can  bo  soon  from  Eq  (78)  that 


C3>i 


1.  _ KV^  -V-  acoslf 


(33) 


Thercforo,  after  some  substitution  Eq  (82)  may  bo 
rev;ritton  in  terms  of  c< as 

- cgAy>VySC.,  ^ 2KVV  - 
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rn 


.(kv% 


^Kcos'^^ 


(8^) 


After  rearranging  and  factoring,  o<  takes  the  final 
form  of 
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Eq  (35)  is  now  used  to  obtain  the  second  derivative  of 
with  respect  to  time.  Initially  this  equation  is 


^ A [c^v ) 
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Rv  + kv  -^kV\cos'K  - cosv  -+- 

<^K  W Sm-if-if  + V (2K  - f~  ®'e'j  + 


(86) 


After  factori.na  tt^ia  equation  can  be  written  as 


o<^  =■  <=><■ 


^ * ** 
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Now  that  the  time  derivative  equations,  bq  (88)  and  Eq 
(37),  are  derived  those  quantities  can  be  used  in  further 
deve  1 opr.ien t of  the  nonlinear  equations.  It  must  be 
recalled  tliat  these  angle  of  attack  derivatives  were 
derived  from  an  estimated  angle  oi  attack  equatiott, 

Eq  (76). 

Tl^e  Moment  Equation 

For  an  aircraft  in  flight  the  moment  equation  is 


(00) 
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Where 


ana 


e = 


I 


Eq  (30)  can  be  rewritten  a; 


JVi 


M 


S c 


(09) 


But  the  pitching  moment  coefficient,  C^,  can  be  expanded 


M 


TT  (s^, 


(90) 


Eq  (90)  can  be  substituted  in  Eq  (09)  to  yield 
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(91) 


This  equation,  Eq  (91),  is  the  second  of  two  equations 
forming  a set  of  simultaneous  equations  which  arc  used 
for  computing  elevator  deflection,  ^ , and  angle  of 
attack,  c>^  . The  first  equation  was  Eq  (79), 


F’eforeriCG  ShahoG 

With  all  the  prGlitninary  equation  derivations 
completed,  the  seven  reference  state  equations  will  be 
developed  in  this  section.  The  seven  reference  states 
comprisinq  the  reference  state  vector  are: 

reference  velocity  in  the  x-direction 

v;,-,  reference  normal  velocity  in  the  z-direction 

I'v 

Pp,  reference  pitch  rate 
reference  pitch  angle 
reference  elevator  deflection 
reference  thrust  (T^  is  equation  notation) 

hj  reference  or  desired  altitude 
d 

Reference  Thrust.  For  computation  of  reference 
thrust  the  forces  are  summed  along  the  velocity  axis  a 
shown  in  Pig.  17.  The  force  equation  is 


Tp,  cos  oc  - D — N\g  s’m 


(9 


Figure  17.  Forces  in  Velocity  Direction 
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Reference  Pitch  Rq^.c.  The  reference  pitch  rate  is 
derived  from  the  basic  definition  of 


Using  the  reference  angles 

« 

Where  oi  r,  comes  from  Eq  (35) 
and  comes  from  Eq  ('IG). 

Reference  Pitch  Annie. 
©,  is  derived  also  from  t'n 


(100) 

(101) 

with  ot „ reolaced  bv  ^ 

The  reference  pitch  angle, 

! basic  definition 


v^hero  <5<  is  from  Eq  (97)  and  is  from  Eq  (43), 

Reference  Altitude,  The  reference  or  desired  alti- 
tude, h^ , is  computed  in  the  reference  path  generator 
and  is  supplied  to  the  nonlinear  reference  state  cenorator- 
Sumrnary . The  seven  reference  r>tatc3  and  their 
respective  equations  arc  listed  below: 

u,  = V COS  o4g 

" V sm  °<R 

S,(  - tq  (96) 


G1 


- Eq  (9-^) 

- reference  path  generator 


Truth  Model 

In  order  to  calculate  the  actual  values  of  each  of 
the  SGV^n  states,  a truth  model  is  required.  This  section 
provides  sets  of  equations  which  calculate  the  actual 
values  for  the  seven  state  variables.  Since  the  M0D6DF 
program  (see  Appendix  D)  integrates  the  soccified  equa- 
tions, various  first-order  differential  equations  are 
dev'^loped . 


Basic  Equations.  The  force  equation,  Eq  (92),  is 
again  used  except  it  is  solved  in  terms  of  V as 

T COS<x  ~ f ^ "*  ^ 


V = 


(103) 


m 


This  equation  for  V can  bo  integrated  to  solve  for  actual 
velocity,  V^. 

The  normal  force  equation  can  be  written  (Ref  17:207) 


which  can  be  solved  for  'if  as 

5 SC  - 


V-.  ““  innaCcs'ilJ  + c< 


( 104) 


V - 


CoS 


f!  +T  sif\  Ot 


( 105) 


hf\ 


Eq  (105)  may  be  integrated  to  yield  the  actual  flight 
path  angle, 
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The  moment  equation  can  be  used  to  solve  for  the 
pitch  rate  derivative.  This  equation,  previously  v;ritton 
in  the  form  of  Eq  (83),  may  bo  rev/ritten  as 

* Q S C. 

(106) 

y 

Eq  (106)  may  be  integrated  to  compute  the  actual  pitch 
rate,  q^.  Then  can  be  used  to  find  the  actual  pitch 
angle  since 


G 


(107) 


Referring  back  to  Chapter  V,  Fig.  12  illustrates 

« « 

the  relationship  betv/oen  h,  R and  V,  Previously  derived 


equations  for  range  and  altitude  derivatives  are 

R — VcOsX  (37) 

and 


U = Ni  Svrv  'i  (36) 

Integration  -f  Eqs  (37)  and  (36)  result  in  actual  range, 
R^,  and  actual  altitude,  h^,  respectively. 

The  actual  elevator  deflection,  is  computed  by 

integrating 


(108) 
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where  7^  is  the  time  constant  and  is  the  reference 

elevator  deflection. 

Likewise,  the  actual  thrust  is  found  by  integrating 


“7^ 

'P 


( 109) 


where  'T'  is  the  time  constant  and  is  the  reference 

thrust , 

Actual  States.  For  convenience,  the  equations  for" 
each  of  the  seven  actual  states  are  summarised  below: 

For  velocity,  t) , the  equation  is 

(UO) 

v/here  is  actual  velocity  and  is  computed  from 

(111) 

Normal  velocity,  w,  is 


= ©A  " 


W.  r:  M.  sin 

A A 


For  pitch  rate,  q,  the  equation  is 


( 1 ,1  P ) 


(113) 


For  pitch  angle,  e , ttio  equation  i: 

% 


G.  = I a 


A 


< 11^1) 


04 


nfnfi  I 


ppf^gmMLiK 


py  " TTT^.r  *• 


The  actual  elevator  deflection  is 


The  actual  thrust  is 


u 


And  the  actual  altitude,  It,  is 


sm  it 


(115) 


(116) 


(117) 


Initial  Conditions.  In  order  to  inteyiate  each  of 
the  first  order  differential  equations  some  initial 
conditions  must  be  set.  These  are  designated  below  for 
each  of  the  quantities  which  is  integrated: 


=0  rad 
o 

V = 647.3  ft/scc 
o 

q =0  rad/soc 

h = 500  ft 
o 


= .066  rad 
- -.065  rad 
= 4729  lbs 
= 0 ft 


Comnutat ion  of  Updated  Control 

fJow  that  the  reference  states  and  the  actual  states 
have  been  computed  it  is  appropriate  to  use  both  to  find 
the  new  control.  With  reference  to  Pig,  13  of  this 
chapter,  the  state  deviation  vector,  can  be  computed 
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by  differencing  the  rofc-ronce  and  truth  state  vectors. 

The  feedback  gain  matrix,  G,  has  been  computed  from  tnc 
linear  model  in  Chapter  III  and  can  be  multiplied  by  the 
state  deviation  vector,  that 

W - G^x  <112) 


Eq  (IIG)  is  the  control  law  which  calculates  the  feedback 
commands  for  the  elevator  and  engine. 

Both  the  total  elevator  control,  C,  and  the  total 
thrust  control,  A , arc  computed  by  addition  of  the 
respective  feedback  command  to  the  reference  command 


and 

A = + ?A 


(119) 


(120) 


This  update  of  the  two  controls  then  completes  one  cycle 
of  the  nonlinear  model  computation  loading  to  a new  set 
of  controls.  The  new  control  values  are  now  used  to 
adjust  the  reference  states  and  the  actual  states  as  the 
simulation  progresses  in  time  along  the  range  direction. 


VII 


Results 


This  chapter  presents  the  results  of  the  design  study 
of  the  proposed  optimal  terrain  follov/ing  controller  using 
a cubic  spline  reference  path.  The  chapter  is  organised 
into  three  main  sections.  The  first  section  shows  the 
resxjltant  numerical  data  for  the  system  matrices.  The 
second  section  illustrates  the  performance  of  the  linear- 
ized model.  And,  in  the  third  section  an  evaluation  is 
made  of  the  nonlinear  controller  model. 

System  r^.atrices 

The  open  loop  dynamics  matrix,  A,  and  the  control 
distribution  matrix,  B,  have  been  discussed  previously. 

In  this  section  the  actual  numerical  values  of  each  matrix 
are  presented.  The  A matrix  of  Eq  (27)  is  shov/n  in  Fig, 
20.  The  3 matrix  of  Eq  (28)  is  illustrated  in  Fig.  19 
below. 


0,  0. 

0.  0. 


10.  0. 

0.  1. 

l_  o-_ 

Figure  19.  The  Control  Distribution  Matrix, 
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The  Output  Matrix 


In  Chapter  III  the  vector 


was  shown  to  be  the  output  deviation  vector.  The  output 
deviation,  Sv,  for  this  thesis  consists  of  the  following 
four  output  variables; 

(a)  Sh  - altitude  deviation 

(b)  ^h  - vertical  velocity  deviation 

(c)  Sh  - vertical  acceleration  deviation 

(d)  $u  ->  forward  speed  deviation 

The  two  variables  Sh  and  ah  v/ere  sized  by  dividing  these 

2 

terms  by  and  respectively.  The  four  output  variable 
equations  are  shown  belov;. 


Sk  (121) 

it 

V.  Vo  \ ‘ 

% ^ (12.) 

H)  \ 


(124) 


These  four  output  equations  can  be  written  in  the  form  of 
Eq  (6)  in  Chapter  III  as 
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0. 

0. 

0. 

0. 

0. 

0. 

0 

-1 

Vo 

0, 

0. 

'iif 

-Zw 

-Zi. 

V/ 

0. 

o 

-2o 

Vo'- 

1. 

0. 

0. 

0. 

0. 

O. 

1 

0. 

Sx 

0. 

(125) 

n 


When  the  numerical  values  are  substituted  the  trans- 
formation matrix,  C,  used  for  this  thesis  is  shovm  in 

r icj  . X , 

Weighting  Matrices  Q and  R.  The  method  for  deter- 
mining the  Q and  R weighting  matrices  was  discussed  in 
Chapter  III,  Doth  matrices  were  developed  in  diagonal 
form.  The  final  forms  of  the  Q and  R matrices  arc  shov/n 
in  Figs.  22  and  23.  These  values  resulted  from  the 
iterations  described  in  the  next  paragraph. 

The  Di'yson  and  llo  method  (Ref  4:119)  was  used  to 
estimate  reasonable  values  for  the  diagonal  elements^ 
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Then  these  values  were  perturbed  3.n  attonpting  to  achieve 
good  performance  for  the  controller.  The  corresponding 
maximum  allowable  deviations  associated  with  the  Q and  R 
diagonal  elements  are  shown  in  Figs.  2*1  and  25. 

State  Deviation  V/eiohtino  Matrix.  In  Chapter  III 

the  cost  function  v;as  'Written  in  terms  of 

This  required  the  state  deviation  weighting  matrix,  P, 

T 

which  was  shown  to  be  equivalent  to  C QC.  This  matrix 
v/ith  its  numerical  values  is  shown  in  Fig.  26. 

Computations  of  CPTSYS.  Using  all  the  system  matrices 
and  the  vjoighting  matrices,  the  OPTSYS  program  calculated 
the  optimal  steady-state  gains  for  the  terrain  following 
controller  and  computed  the  eigenvalues  of  the  system. 

The  steady-state  Ricatti  gain  matrix,  2,  is  presented 
in  Fig.  27  and  the  feedback  gain  matrix,  G,  is  shovni  in 
Fig.  20. 

The  seven  eigenvalues  of  the  open  loop  dynamics 
matrix.  A,  are  shown  in  Fig.  29.  All  the  open  loop 
eigenvalues  are  stable. 

The  closed  loop  dynamics  matrix,  A+DG,  results  in 
eigenvalues  which  have  been  moved  further  into  the  stable 
region  of  the  complex  plane.  This  increase  in  stability 
is  shov;n  by  comparison  with  the  corresponding  open  loop 
eigenvalues  in  Fig.  29.  The  dominant  complex  roots  pro- 
duce a damping  factor  close  to  .707.  This  dampirig  factor 
was  desired  so  that  the  system  would  possess  good  time 
response  and  damping  v/ith  minimal  peak  overshoot.  It  was 
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also  required  that  tlie  system  not  be  overly  damped  since 
that  v^ould  cause  excessive  acccleraticn  and  undue  crow 
discomfort. 


Stability  of  the  Lincaricod  Perturbation  dauations 

A computer  program  v/as  developed  to  verify  the 
stability  of  the  linearized  perturbation  equations  derived 
in  Chapter  17.  Since  the  coefficients  in  these  equations 
determine  the  value  of  the  A matrix  and  ultimately  dictate 
the  optim^.l  feedback  gains,  it  was  important  to  insure 
their  stability.  The  matrices  A,  3,  G,  and  A+BG  v;cre 
used  to  simulate  the  system  in  closed  loop  form. 

Since  the  linearised  equations  are  in  perturbation 
form  a straight  and  level  flight  profile  was  assumed  as 
the  zero  reference  condition.  Initial  conditions  were 
set  on  one  of  the  seven  perturbation  stats  variables 
v;hile  the  remaining  six  were  set  to  zero.  Then  the  seven 
differential  equations  (see  Eqs  (19)  through  (25;)  were 
integrated  from  0 to  20  seconds.  In  addition,  Eq  (122) 
was  inteqrated  to  check  the  output  vertical  velocity 
performance.  Figures  30,  31,  32,  and  33  illustrate  the 
linoarix.cd  equations  rjorformance  for  an  initial  perturba- 
tion iri  x-dircction  velocity  of  25  ft/scc.  Figures  3f  to 
3?  correspond  to  o.n  initial  picrhurbation  in  altitude  of 
25  feet.  And,  Figures  38  to  show  how  a perturbation 
of  200  lbs  of  thrust  affects  the  closed  loop  linear  model. 

Also  included  in  Figs.  42  to  45  arc  plots  of  t)ic 


E'V'h 

|g  'jt;, 


VERT  VEL  OEV  ALT  DEV  (FT) 

.40  0.00  0,40  0.80  0-0'7  0.15  0.23 
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Figure  40,  Rcr.ponrc  to  Thrurt  Pcrtur’oat  i on  of  200  lbs 


stoady-sLate  pei'f orniance  of  the  controller  responding  to 
ai'i  initial  poi  turbat  ion  ot  25  ft/soc  in  x-direction 
velocity.  These  plots  show  that  in  steady-stat  e t tie 
controller  restores  all  states  to  the  oero  flight 
cond  it  ions . 

In  all  cases,  tlic  time  response  or  ttie  liricar  system 
is  shown  to  be  qnite  good.  There  arc  no  large  fluctua- 
tions in  response  and  ttie  system  is  restored  to  the  cero 
reference  in  a reasonable  amount  of  time. 

Performance  of  the  Nonlinear  Con  troll er 

Much  of  tlio  initial  study  of  this  tlicsis  consisted 
of  verification  of  the  equations  used  in  the  nonlinear 
reference  and  truth  models.  Once  these  equations  and  the 
associated  computer  programs  were  investigated  it  was 
possible  to  evaluate  the  performance  of  the  controller 
with  the  nonlinear  reference  and  truth  models. 

The  M0D6DF  program  was  vised  to  compute  a simulated 
reference  path  for  the  controller.  This  simulated  path 
consisted  of  a simpile  pvillvip  maneuver  in  which  an  air- 
craft climbs  over  an  obstacle  and  levels  off.  This  is 
shown  in  Fig.  TG.  For  this  simple  pullup  maneuver  th.e 
acceleration  and  g rate  limits  can  be  varied  and  the 
duration  of  the  pullup  can  be  specified. 

Tv;o  different  reference  generation  schemes  were 
developed.  One  involved  a cubic  spline  path  and  the 
other  was  a quartic  (polynomial  of  degree  four)  spline 
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path.  Hodulcb  wore  designed  tor  program  MOD6UI'  which 
consisted  of  a spline  path  generation  routine  and  a 
corresponding  iiitfu'polaL  iot'i  routii'^e.  These  reference 
path  generation  sudroutines  were  then  linked  to  ttie 
nonlinear  controller  modc‘1  to  evaluate  pei'i  ot  mancc. 

For  the  ijuartic  Sf^line  rofci'enco  ['ullup  nanouver  ttie 
initial  altitude  was  500  feet.  Over  a range  of  10,000 
feet  tt'.e  aircraft  climt-ed  1090  leet  to  level  off.  The 
spo'cified  limits  on  k (see  Fig.  -16)  were  .5g  for  this 
profile.  Ttio  resulting  quartic  reference  path  was  tracked 
by  tlic  controller  wi.th  a maximum  altitude  error  of  G fcot- 
Thu  velocity  deviation  vs'ar.  loss  tliat^  2 ft/scc:  at  level 
off.  Both  the  altitude  and  velocity  tracking  errors  were 
regulated  to  zero  after  completion  of  ttie  pullup  maneuver. 
Tlio  velocity  error  settling  time  was  slightly  longer  than 
the  altitude  settling  time  after  1-evel  off. 

The  cubic  spline  pullup  maneuver  was  also  performed 


for  a 

.5g  limit 

on  k . Over 

a range 

of 

9,000  feet 

the 

cubic 

spline  lei 

Ci  (.fi  po  t I'l 

1.  n v,,  I o 

i Q 1 1 i.  ii  1 

ijijo 

f L'om 

500  f 

cet  to  1280 

feet.  The 

opt imal 

contc  oiler 

tracked 

the  reference  altitude  within  4 . ‘3  feet.  The  velocity 
tracking  error  for  ^he  cubic  sp^line  reference  pullup  was 
less  than  2 ft/sec.  Like  the  quartic  pullup  maneuver, 
both  the  velocity  error  and  the  altitude  error  were  nulled 
after  level  off  by  the  feedback  controller. 

The  optimal  controller  was  next  flown  using  the 
optimized  cubic  spline  reference  path  described  in 
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Chapter  V,  Over  the  first  15,000  feet  of  this  terrain 
profile  the  standard  deviation  in  altitude  was  160  feet. 
The  optimal  controller  was  capc'.blc  of  maintaining  the 
altitude  tracking  error  at  less  than  1,5  feet  and  kept 
the  velocity  tracking  error  less  than  1,5  ft/sec.  This 
was  very  good  performance  even  though  the  terrain  was 
rather  benign.  The  thrust  fluctuation  was  reduced  when 
compared  to  the  Boeing  optimal  controller  flying  over  a 
similar  low  level  terrain  profile.  But  the  actual 
thrust  did  lag  the  required  reference  thi'ust  because  of 
the  simplified  thrust  control  systera  developed  for  this 
thesis. 

In  all  the  flight  profiles  flown  by  the  optimal 
terrain  following  controller  so  far,  the  tracking  errors 
have  not  been  large  and  can  be  considered  quite  accept- 
able. Wnenever  an  altitude  tracking  error  exists  the 
actual  aircraft  altitude  has  always  been  greater  than 
the  reference  altitude  thus  promoting  a safe  terrain  fol- 
lowirig  mission.  Even  though  these  test  profiles  have  not 
been  loiig  in  duration  they  do  show  the  potential  for  good 
controller  performance  using  an  optimized  cubic  spline 
reference  path. 
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VIII. 


Cone  ] un  ion.'i  arui  :^ecoinnendo  1 i on  5^ 


The  feedback  controller  ic  a stable  controller 
capable  of  track  inn  a cubic  spline  reference  path.  The 
cubic  spline  generator  produces  a reference  patli  possessing 
good  tracking  quality.  The  resulting  controller  provides 
good  path  control  and  only  slight  tracking  errors  in  alti- 
tude and  velocity  provided  the  state  deviations  remain  in 
the  linearised  region  aaout  the  nominal  reference  trajec- 
tory. By  allowing  variations  in  velocity,  the  thrust 
fluctuation  is  reduced.  This  is  essential  to  increasing 
engine  life  and  decreasing  fuel  consumption. 

The  v.'oighting  matrices,  Q and  R,  are  very  influential 
in  the  performance  of  the  terrain  following  controller. 
These  matrices  can  be  judiciously  adjusted  to  strengthen 
the  coi'iti.'oller  operation  and  cause  the  aircraft  to  fly 
closer  to  the  reference  path. 

The  following  recommendations  result  from  v/ork  on 
this  design  study 

(a)  investigate,  ttirough  more  sophist icated  tech- 
niques like  root  square  locus,  the  eigenvalue 
placement  due  to  varying  the  Q and  1^  matrices 

(b)  evaluate  this  nonlinear  controller  over  a 
greater  variety  of  terrain  profiles 

(c)  investigate  ttic  linearising  assumptions  made  in 
Chapter  IV  and  determine  v;hctlier  a lotating 
coordinate  reference  frame  tangent  to  the 
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reference  trajectory  is  mere  accurate 

(d)  develop  a nore  sophisticated  thrust  and  velocity 
control  system  for  improvinq  the  controller 
performance 

(e)  modify  the  constant  energy  assumption  to  allow 
variations  in  energy  within  a specified  energy 
envelope 


'A 
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Appendix  A 


Constants  of  the  Aircraft  Model 


This  appendix  contains  all  the  modelling  constants 
used  in  the  controller  design  simulation.  The  data  rep- 
resents a large  fighter-bomber  aircraft  flying  at  sea 
level.  The  conditions  for  a standard  atmospheric  day  are 
assumed. 


Density 

Speed  of  sound 
Nominal  aircraft  velocity 
Aircraft  weight 
Wing  area 


^ = 2.378E-03  slug/ff^ 

A = 1116  ft/sec  (Ref  17:465) 


V = 647.3  -Ft/sec  (Mach  .58) 
o 


w 


55000  lbs 
2 


S « 647  ft 


Wing  mean  aerodynamic  chord  c - 18,678  ft 
Dynamic  pressure  at  nominal  velocity 

q = S = 322307.6  slug-ft/sec^ 

Moment  of  inertia  about  aircraft  y-axis 
I = 380000  slug-ft/sec^ 

y 
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Appendix  B 


Stability  Derivatives 

This  appendix  presents  all  the  stability  derivatives 
in  the  controller  design  and  explains  their  derivation. 
First  of  all,  the  dimensional  stability  derivatives  are 
designated.  Then  the  non-dimensional  stability  deriva- 
tives are  explained.  And,  finally,  the  numerical  values 
of  all  the  stability  derivatives  are  compiled. 

Dimensional  Stability  Derivatives 

The  linear  model  in  Chapter  IV  illustrated  how  the 
linearized  equations  of  motion  contain  dimensional 
stability  derivatives.  These  dimensional  stability  de- 
rivatives are  the  coefficients  in  the  open  loop  dynamics 
matrix,  Eq  (271,  These  dimensional  derivative  values 
result  from  data  computed  by  ASD/XRHA,  Wr ight-Patterson 
AFB,  Ohio.  The  definition  of  each  dimensional  derivative 
is  included  with  its  corresponding  symbol.  In  addition, 
for  dimensional  derivatives  required  for  further  com- 
putation, the  corresponding  equation  and  units  are 
included.  The  dimensional  derivative  equations  and  the 
units  are  compiled  from  Reference  16. 


Dimensional  l/ariation  of  x-force  with  Speed 


VVT,  \L 


{ sec  ) 


(B-1) 


105 


, Dimonsional  Variation  of  x-force  with  Elevator 
Dcf lec  t ion 


(ft/sec^) 


(B-2) 


X. , Dimensional  Variation  of  x-force  with  Anqle  of  Attack 


- ^ ^ (ft/sec^) 


(B-3) 


lv\ 


Z.,,  Dimensional  Variation  of  c-force  with  Speed 


2.- 


- a(r 


^ - Wo  / 


(sec“'^) 


(Q-4) 


)rr\ 


u 


Z^,  Dimensional  Variation  of  z.-force  with  Elevator 
Del  lection 


7 - 

H' 


S C|_^ 

A b 


{ f t/sec^) 


V a-3  j 


Z^,  Dimensional  Variation  of  s-force  with  Anqle  of  Attack 


Z.= 


ft/SGC  ) 


( 3-6 ) 


hrx 


Dimensional  Variation  of  Pitching  Moment  with  Speed 
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( f t-sec) 


(B-7) 


~1 


, Dimensional  Variation  of  Pitching  Moment  with  Elevator 
Def lection 


1 ^ ^ 


I 


1 


( sec~^) 


(B-8) 


Dimensional  Variation  of  Pitching  Moment  with  Angle 
of  Attack 


^ ScC 


I 


( sec"*^) 


(B-9) 


The  nine  equations,  Eq(B-l)  through  Bq  (6-9),  are 
required  to  compute  the  non-dimensional  stability  deriv- 
atives shown  later  in  Appendix  B. 

It  can  be  seen  that  the  open  loop  dynamics  matrix, 

A,  in  Eq  (27)  does  not  include  the  dimensional  derivative 
terms  X^,  Z^,  or  coefficients.  Instead,  the  longi- 

tudinal equations  of  motion  include  X , Z , and  M 
respectively.  With  a simplified  assumption  these  three 
varianles  in  w can  be  transformed  into  variables  in  «»<■ 
for  computation  purposes.  As  seen  in  Fig.  B-1  for  small 
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Figure  B-1.  Velocity  Relationships 


angles  of  attack, o<,  the  following  equation  is  true 
CRef  12:72,  6:11) 


^ V. 

There  is  also  a relationship  between 
derivatives  as  follows 


dimensional 


tB-lO) 


stability 


(B-11) 


Substitution  of  Eq  (3-10)  into  Eq  (B-11)  results  in  an 
important  relation 


X 


(B-12) 


Now,  since  the  value  of  is  a known  quantity  from  the 
open  loop  dynamics  matrix,  it  is  possible  to  solve  for 


8 


using  Eq  (B-IZ).  There  Is  an  available  equation,  Eq  (B-3), 
which  is  written  in  terms  of  so  it  car  be  used  to 
calculate  the  non-dimensional  stability  derivatives.  The 
same  process  developed  for  is  used  for  the  and  the 
equations.  This  is  the  method  needed  to  solve  for 
several  unknowns  in  the  nonlinear  model  equations., 

Non-Dimensional  Stabili ty  Derivatives 

The  non-dimensional  stability  derivatives  are  not 
used  in  the  open  loop  dynamics  matrix,  A,  but  are  required 
for  the  nonlinear  model  equations.  The  non-dimensional 
stability  derivatives,  with  their  definitions  and  perti- 
nent equations,  are  listed  below  (Ref  16); 


is  the  lift  coefficient 


c 


L 


C.  is  the  aircraft  lift  curve 


\r 


oC 


C,  is  the  control  surface  lift  effectiveness 


IS 


(B-13) 


(B-14) 


(B-15) 
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C,  is  the  variation  in  lift  coefficient  v/ith  forward 
- U 

speed 


“-u  - 


(B-16) 


is  the  lift  coefficient  for  zero  angle  of  attack. 
Cjj  is  the  drag  coefficient 


1 


(B-17) 


C_.  is  the  variation  in  drag  coefficient  with  angle  of 
attack 


’0  . ^ 


\ r 


(B-18) 


oL 


is  the  variation  of  drag  coefficient  with  elevator 
deflection 


C- 

“s 


U 


Cp^  is  the  variation  of  drag  coefficient  with  forward 
speed 


>C, 


■*>«  - 


(B-20) 
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is  the  drag  coefficient  for  .-ero  angle  of  attack, 
o 

C is  the  pitching  moment  coefficient 
m K -a 


c 


(B-21) 


C is  the  variation 
angle  of  attack 


C is  the  variation 
elevator  deflection 


in  pitching  momeiit 

of  pitching  moment 


coefficient  with 

(B-22) 

coefficient  with 


c 


iC 


(B-23) 


C is  the  variation  in  pitching  moment  coefficient  with 


my 

forward  speed 

r 


yrr\ 


- 


KB-ZA) 


C is  the  variation  in  pitching  moment  coefficient  with 
m 

q 

pitch  rate 


) 


(B-25) 


C is  the  pitching  moment  coefficient  for  zero  angle  of 
attack. 


111 


All  of  these  non-dimensional  stability  derivatives 
are  required  for  the  nonlinear  model  equations.  The 
numerical  values  are  shown  in  the  next  section. 

Stability  Derivative  Data 

The  values  of  the  dimensional  stability  derivatives 
are  known  from  the  open  loop  dynamics  matrix,  A (see 
Chapter  IV).  From  these  known  values  and  the  model 
constants  all  the  necessary  derivatives  are  calculated. 
The  actual  dimensional  stability  derivatives  are: 


= 

-.007578 

sec"*^ 

.01276 

-1 

sec 

X 

q 

= 

.034 

ft/sec 

= 

-.000314 

f t/sec'^ 

= 

.000713 

sec**^ 

=: 

-.1029 

sec”^ 

= 

-.8198 

-1 

sec 

z 

q 

z= 

-2.175 

f t/sec 

H 

= 

. 

CO 

1 

ft/ sec^ 

= 

-.0000361 

-1 

sec 

-.00013 

(ft/sec ) 

= 

-.006086 

(f t/sec) 

M 

q 

= 

-.201 

-■1 
sec  ~ 

= 

-4.477 

sec-2 

Befo.e  the  non-dimensional  stability  derivatives  are 
listed,  the  zero  angle  of  attack  conditions  must  be  set. 
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These  values  are; 


Cp  - .0106 

o 


With  these  zero  angle  of  attack  quantities,  the  remainder 
of  the  stability  derivatives  may  be  calculated  from  the 
equations  of  the  first  two  sections  of  this  Appenfiix, 


<=L  “ 

.1706 

C,  - .2567 

V 

2,804 

C,  = . 35328 

•',0412 

C_  e .00402 
’ u 

c„  = 

1.6654x10”^ 

C = -.2484 

c « 

-.0053 

C »=  —.2326 

c » 

-.8794 

These  values  are  used 

In  the  variety  of  equations  in 

the 

nonl Ineai 

control ler 

model.  This  data  is  consistent 

with 

actual  flight  test  data  (Ref  7 and  17). 
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Appendix  C 

Calculation  of  Trim  Conditions 

DuiHrig  the  trim  condition  for  stability  the  angular 
rates  arc  zero,  the  linear  rates  arc  v^oro  and,  essentially, 
lift  equals  vjeiqht  (Ref  7:30), 

For  the  trim  condition 


L - i^SC,  = u 


(C»l) 


Substitution  of  the  known  valuer.  Into  Eq  (C-l)  provides 
5 value  for  of 


L. 


i 'irKt 
.1  /VO 


The  corresponding  trim  equation  Is 


ot 


(C-2) 


Eq  (C-2)  is  one  equation  needed  to  solve  for  tlie  two 
trim  angles,*^  and  S. 

The  pitching  moment  coefficient  cqi^ation  is  written  as 


^ rw  oC  -f  Ck,  . i 


(C-3) 


For  the  trim  condition  C c 0 and  Fq  (C-3)  bccomci 

m ^ 


o - + c:,„  .c<  + 


(C-4) 


1 M 


From  Appendix  B the  value  of  is  zero  for  the  control- 

o 

ler  design,  so  Eq  (C-4)  reduces  to 

0 “ ^ (C-5) 

This  is  the  second  equation  required  to  solve  for  the 
tv;o  unknown  trim  angles.  Eq  (C-2)  and  Eq  (C-5)  can  be 
solved  simultaneously  to  yield  the  trim  conditions  of 

raJ 

and 

^ = -.0582  nJ 

TRII1 

These  two  trim  angles  are  used  for  calculations  in  the 
nonlinear  reference  model. 
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Appendix  D 


MQD6DF  Computer  Program 

The  Litton  M0D6DF  algorithm  is  a s ix-degree-of- 
freedom  digital  simulation  program  useful  as  a tool  to 
analyze  various  guidance  and  control  systems.  It  uses 
a modular  approach  to  FORTRAN  programming  which  allows 
the  user  flexibility  to  adapt  modules  to  various  physical 
systems*  This  appendix  consists  of  a summary  of  th^ 
M0D6DF  program.  For  more  detail  about  the  algorithm 
the  reader  should  consult  Reference  5. 

Program  Structure 

The  M0D6DF  program  is  structured  for  systems  engi- 
neers rather  than  programmers.  It  permits  an  engineer 
to  use  the  overall  program  without  being  too  concerned 
about  numerous  programming  details.  An  engineer  only 
needs  to  modify  an  independent  module  representing  the 
physical  system  of  interest. 

The  basic  structure  of  the  M0D6DF  program  is  built 
around  three  main  levels  of  subroutines  (see  Figure  D-l) . 
These  three  types  of  subroutines  are: 

(a)  Type  0 - Executive  subroutines  which  manage 
and  supervise  the  overall  MODGDF  program 

(b)  Type  1 - Operational  subroutines  which  pr-ovide 
control  at  the  beginning  of  each  integration 
step 
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(c)  Type  2 - Modules  which  are  constructed  by  the 
user  according  to  the  physical  system  being 
invest igated 

The  executive  and  operational  subroutines  are  not 
normally  modified  by  the  user.  The  modules  are  the 
"black  boxes"  of  the  M0D6DF  program  which  represent 
modelling  equations  from  the  physical  system.  The  proper 
development  of  these  modules  is  the  main  concern  of  an 
engineer  utilizing  the  MOD6DF  program. 

Operat ion 

The  M0D6DF  program  operates  basically  as  an  inte- 
gration subroutine  package.  It  utilizes  a fourth-order 
Runge  Kutta  integration  scheme  in  closed  form,  Vjithin 
each  module  there  is  a number  of  differential  equations 
represent irg  the  physical  system.  The  user  identifies 
the  locations  of  the  variables  which  are  integrated 
by  v/riting  certain  specifications  in  an  Initialization 
module.  Uriug  the  integration  alcjorithm,  the  system 
solutio''  progresses  at  each  point  in  time  by  processing 
the  specified  d: i f erential  equations  to  advance  the  solu- 
tion one  step  ^t  in  time.  Computed  values  of  various 
variables  mejy  be  stored  at  each  stop  At  for  later  use 
as  output  data  or  for  plotting.  The  program  can  also 
generate  raiidom  noise  for  Monte  Carlo  simulation  and  can 
perform  covariance  simulation. 

Most  o"  the  communication  betv;een  subroutines  is 
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accomplished  through  COMMON  locations.  Many  of  the 
COMMON  storage  locations  are  specified  (Ref  5:23),  but 
others  are  reserved  for  the  variables  in  the  system 
modules.  Therefore,  the  user  docs  not  have  to  use  CALL 
statements  or  argument  lists  when  performing  analysis. 

The  COMMON  locations  are  arranged  so  that  each  variable 
being  integrated  is  associated  with  four  consecutive 
COMMON  locations.  The  first  location  is  for  the  deriva- 
tive of  the  variable  which  is  specified  by  the  differential 
equation  in  the  system  module.  The  second  and  third 
locations  v;ere  previously  for  the  lower  and  upper  bounds, 
respectively,  on  the  integration  error  of  that  variable; 
butc  these  are  no  longer  in  use.  And,  the  fourth  loca- 
tiofi  is  for  the  variable  itself.  Figure  D-2  illustrates 
an  example  of  how  this  COMMON  arrangement  works  for  the 

Variable 

V 

V 

Figure  D-2.  COMMON  locations  for  a typical  integration 
basic  equation  F = mV 

Input  Cards 

The  input  cards  to  MOdGDF  are  of  ten  types.  These 
input  cards  are  closely  related  to  the  operational  sub- 
routines and  the  system  modules.  As  a result,  the  input 
cards  essentially  dictate  how  the  MODGDF  program  will  be 


Value 

F/m 

^F/m 


COt-lMON  location 
301 
304 


rj 


"1 


» 


'I 


I 


'1 

i 


i 
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implemented  for  any  particular  data  run.  There  is  a 
basic  format  for  the  input  cards  which  must  be  strictly 


followed.  The  ten  types  of 

Tx£e 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 


input  cards  are  shown  helow. 

Input  Description 

Operational  subroutines  to 
be  used 

System  modules  to  be  used 

Numerical  input  data  includ- 
ing initial  conditions 

Variables  to  be  printed  out 

Variables  to  be  controlled 

Parameters  of  random  number 
generators 

Variables  to  be  plotted 

Spare 

Spare 

Indicates  end  of  a data  run 


r ^ 
I 


-3 

-3 
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Appendix  E 


The  OPT SYS  Proqram 


This  appendix  explains  the  input  and  output  para- 
meLers  of  the  OPTSYS  computer  program.  It  does  not 
describe  \;he  program  operation  in  detail  since  that  is 
the  purpose  of  Reference  3, 


Program  Description 

OPTSYS  is  a digital  algorithm  which  utilizes  optimal 
control  theory  to  calculate  optimal  gains  for  a system  of 
the  general  form 


(E-l) 


using  a cost  function  v/ritten  as 


J = if  (s/PSx  + SiJRWyi 


(E-2) 


The  OPTSYS  program  sol''-es  the  steady-state  form  of 
the  matrix  Ricatti  equation  to  obtain  the  corresponding 
Picatti  gain  matrix,  S.  It  is  possible  to  compute  the 
feedback  gain  matrix,  G,  using 


G c -r’b"5 


(E-3) 


where  all  three  matrices  on  the  right  hand  side  of 
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Eq  (E-3)  are  known.  The  linear  control  law 


which  permits  substitution  of  Eq  (E-4)  into  Eq  (E-1)  to 
yield 

Sx  = (A+B&')  Sx  ((..3, 

The  resulting  matrix,  A+BG,  is  called  the  closed  loop 
dynamics  matrix.  The  stabilicy  characteristics  of  this 
closed  loop  matrix  determines  how  well  a feedback  con- 
troller operates. 

Input  Parameters 

OPTSYS  requires  the  following  parameters  for  opera- 
tion ; 

(a)  A the  open  loop  dynamics  matrix 

(b)  B the  control  distribution  matrix 

(c)  P the  state  weighting  matrix 

(d)  R the  control  weighting  matrix 

Using  these  four  matrices  which  correspond  to  the  standard 
form  of  Eqs  (E-1)  and  (E-2),  the  optimal  gains  are  calcu- 
lated. 

Output  Parameters 

The  output  quantities  of  the  OPTSYS  program  are: 
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(a)  the  eigenvalues  and  eigenvectors  of  the  open 
loop  dynamics  matrix,  A 

(b)  t^le  steady-state  Ricatti  gain  matrix,  S 

(c)  the  feedback  gain  matrix,  G 

(d)  the  closed  loop  dynamics  matrix,  A+DC,  and  its 
associated  eigenvalues  and  eigenvectors 

Additional  Uses 

For  the  ourpose  of  this  thesis  OPTSYS  performed 
calculations  as  a deterministic  controller.  It  is 
possiole  to  introduce  noise  in  OPTSYS  and  to  perform 
stochastic  calculations  uijng  estimation  theory.  Proce- 
dures for  this  mode  of  operation  may  te  found  in 
Referen.r:e 
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